#=============================================================================#
#
# PROJECT:        Who Pays for Peace?
# AUTHORS:        ** anonymized for review **
# CONTACT:        ** anonymized for review **
# LAST MODIFIED:  July 5, 2022
# 
#=============================================================================#
#
# This R file contains the code used to replicate the average treatment effects 
# (at Wave 1)
# 
#=============================================================================#


# Initial settings ------------------------------------------------------------

#rm(list=ls())
getwd()

## Install and load all necessary packages ------------------------------------
# ipak function: install and load multiple R packages.
# check to see if packages are installed. Install them if they are not, then load them into the R session.

ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("tidyverse", "ggplot2", "ggpubr", "jtools", "effectsize") 
ipak(packages)


## Load the data --------------------------------------------------------------
mydata <- readRDS("data/T1-clean-data.rds")


# Means per condition ---------------------------------------------------------

#means, sd, and CI per condition, compensations
sum_c <- mydata %>% 
  group_by(comp_condition) %>%
  rename(condition = comp_condition)  %>%
  summarise(mean = mean(compensation, na.rm = TRUE),
            sd = sd(compensation, na.rm = TRUE),
            n = n()) %>%
  mutate(se = sd / sqrt(n),
         lower.ci = mean - qt(1 - (0.05 / 2), n - 1) * se,
         upper.ci = mean + qt(1 - (0.05 / 2), n - 1) * se,
         outcome = "Compensation")

#means, sd, and CI per condition, war crime punishments
sum_p <- mydata %>%
  group_by(punish_condition) %>%
  rename(condition = punish_condition)  %>%
  summarise(mean = mean(punishment, na.rm = TRUE),
            sd = sd(punishment, na.rm = TRUE),
            n = n()) %>%
  mutate(se = sd / sqrt(n),
         lower.ci = mean - qt(1 - (0.05 / 2), n - 1) * se,
         upper.ci = mean + qt(1 - (0.05 / 2), n - 1) * se,
         outcome = "Punishment")

#means, sd, and CI per condition, both outcomes
sum <- rbind(sum_c, sum_p)
sum


## Figure 1 -------------------------------------------------------------------

#plot preparations
sum$label[sum$condition=="Armenia"] <- "Armenia\n(outgroup)"#nicer labels
sum$label[sum$condition=="Azerbaijan"] <- "Azerbaijan\n(ingroup)"
sum$label[sum$condition=="Both"] <- "Armenia +\nAzerbaijan"
sum$label[sum$condition=="Int. Comm."] <- "International\nCommunity"
sum$label[sum$condition=="Control"] <- "Control\nGroup"
sum$label <- factor(sum$label, levels = c("Armenia\n(outgroup)", "Azerbaijan\n(ingroup)", "Armenia +\nAzerbaijan",
                                          "International\nCommunity", "Control\nGroup"))
pd <- position_dodge(.5) #move points/bars .05 to the left and right
cols <- c("#440154", "#3b528b", "#21918c", "#5ec962", "#fdd525") #colors treatment conditions

#open plot environment
jpeg("figures/fig1.jpg", width = 8, height = 5, units='in', res = 300)

#plot
ggplot(sum, aes(x = label, y = mean, colour = label, shape = outcome)) + 
  geom_errorbar(aes(ymin = lower.ci, ymax = upper.ci), width = .1, position = pd, size = 0.7) +
  geom_point(size = 3, position = pd) + 
  scale_y_continuous(breaks = round(seq(2, 5, by = 0.5),1)) +
  scale_color_manual(values = cols, guide = "none") +
  labs(x="\nExperimental Condition", y="Mean Level of Support\n", 
       shape = "Peace Provision") +
  theme_bw() +
  theme(axis.ticks.length = unit(.25, "cm"),
        axis.text.x = element_text(colour = cols, size = 12),
        axis.text.y = element_text(colour = "gray20", size = 12),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 13))

#save plot
dev.off()


# Average Treatment Effects ---------------------------------------------------

## Compensations --------------------------------------------------------------
lm_c <- lm(compensation ~ comp_condition, data = mydata)#simple OLS regression
summary(lm_c)#check results (unstandardized)
effectsize(lm_c) #standardized coefficients (factors are not standardized!)


## Punishments ----------------------------------------------------------------
lm_p <- lm(punishment ~ punish_condition, data = mydata)
summary(lm_p)
effectsize(lm_p) 

## Plot ATEs
# plot preparations
names(lm_c$coefficients)[1] <- "(intercept)"
names(lm_p$coefficients)[1] <- "(intercept)"
names(lm_c$coefficients)[2] <- "Armenia"
names(lm_p$coefficients)[2] <- "Armenia"
names(lm_c$coefficients)[3] <- "Azerbaijan"
names(lm_p$coefficients)[3] <- "Azerbaijan"
names(lm_c$coefficients)[4] <- "Both"
names(lm_p$coefficients)[4] <- "Both"
names(lm_c$coefficients)[5] <- "IC"

#open plot environment
jpeg("figures/fig-presentation.jpg", width = 8, height = 5, units='in', res = 300)

# plot
plot_summs(lm_c, lm_p,
           omit.coefs = "(intercept)",
           model.names = c("Compensation", "Punishment"),
           colors = c("gray40", "black")) +
  labs(y=NULL, x="\nAverage Treatment Effect (ATE)", 
       shape = "Peace Provision") + 
  theme_bw() +
  theme(axis.ticks.length = unit(.25, "cm"),  
        axis.text.x = element_text(colour = "gray20", size = 12),
        axis.text.y = element_text(colour = "gray20", size = 12),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 13))

#save plot
dev.off()


## Benchmark effect sizes -----------------------------------------------------
# i.e., test role of exposure on support in control group only
control <- mydata[mydata$comp_condition=="Control",]

### Compensations -------------------------------------------------------------
lm_c_exposure <- lm(compensation ~ exposure_2020, 
                    data = control)
summary(lm_c_exposure)
effectsize(lm_c_exposure)

### Punishments ---------------------------------------------------------------
lm_p_exposure <- lm(punishment ~ exposure_2020, 
                    data = control)
summary(lm_p_exposure)
effectsize(lm_p_exposure) 

## Intergroup attribution bias: descriptive evidence --------------------------
prop.table(table(as.factor(mydata$aggression)))
prop.table(table(as.factor(mydata$violence)))
